module eddy_diff_cam

  use shr_kind_mod    , only: i4 => shr_kind_i4, r8 => shr_kind_r8
  use ppgrid          , only: pcols, pver, pverp
  use cam_logfile     , only: iulog
  use cam_abortutils  , only: endrun
  use physconst       , only: gravit, cpair, rair, zvir, latvap, latice, karman
  use diffusion_solver, only: vdiff_selector
  use eddy_diff       , only: ncvmax
  use time_manager    , only: is_first_step
  use physics_buffer  , only: physics_buffer_desc
  use spmd_utils      , only: masterproc
  use phys_control    , only: phys_getopts

  implicit none

  private

  public eddy_diff_readnl
  public eddy_diff_register
  public eddy_diff_init
  public eddy_diff_final
  public eddy_diff_tend

  ! Is UNICON switched on (and thus interacting with eddy_diff via pbuf)?
  logical unicon_is_on

  ! Number of iterations for solution
  integer, parameter :: nturb = 5

  ! Logical switches for moist mixing ratio diffusion
  type(vdiff_selector) fieldlist_wet
  ! Logical switches for molecular diffusion
  ! (Molecular diffusion is not done here.)
  type(vdiff_selector) fieldlist_molec

  integer ntop_eddy, nbot_eddy

  ! Cloud mass constituent indices
  integer ixcldliq, ixcldice

  ! input pbuf field indices
  integer :: qrl_idx   = -1
  integer :: wsedl_idx = -1

  ! output pbuf field indices for UNICON
  integer :: bprod_idx    = -1
  integer :: ipbl_idx     = -1
  integer :: kpblh_idx    = -1
  integer :: wstarPBL_idx = -1
  integer :: tkes_idx     = -1
  integer :: went_idx     = -1

  ! Mixing lengths squared.
  ! Used for computing free air diffusivity.
  real(r8), allocatable :: ml2(:)

  ! Various namelist options to limit or tweak the effects of eddy diffusion.

  ! Pressure defining the bottom of the upper atmosphere for kvh scaling (Pa)
  real(r8) :: kv_top_pressure = 0._r8
  ! Eddy diffusivity scale factor for upper atmosphere
  real(r8) :: kv_top_scale = 1._r8
  ! Eddy diffusivity scale factor for the free troposphere
  real(r8) :: kv_freetrop_scale = 1._r8

  ! The following all have to be set in all cases.
  real(r8), parameter  :: unset_r8 = huge(1._r8)
  ! Maximum master length for diag_TKE
  real(r8) :: eddy_lbulk_max = unset_r8
  ! Maximum dissipation length for diag_TKE
  real(r8) :: eddy_leng_max = unset_r8
  ! Bottom pressure level (hPa) for eddy_leng_max
  real(r8) :: eddy_max_bot_pressure = unset_r8
  ! Moist entrainment enhancement param
  real(r8) :: eddy_moist_entrain_a2l = unset_r8

contains

  subroutine eddy_diff_readnl(nlfile)

    use namelist_utils, only: find_group_name
    use units         , only: getunit, freeunit
    use spmd_utils    , only: masterprocid, mpi_real8, mpicom
    use shr_log_mod   , only: errMsg => shr_log_errMsg

    character(*), intent(in) :: nlfile

    integer unitn, ierr
    character(*), parameter :: subname = 'eddy_diff_readnl'

    namelist /eddy_diff_nl/ kv_top_pressure, kv_top_scale, &
         kv_freetrop_scale, eddy_lbulk_max, eddy_leng_max, &
         eddy_max_bot_pressure, eddy_moist_entrain_a2l

    if (masterproc) then
      unitn = getunit()
      open(unitn, file=trim(nlfile), status='old')
      call find_group_name(unitn, 'eddy_diff_nl', status=ierr)
      if (ierr == 0) then
        read(unitn, eddy_diff_nl, iostat=ierr)
      end if
      if (ierr /= 0) then
        call endrun(subname // ':: ERROR reading namelist')
      end if
      close(unitn)
      call freeunit(unitn)
    end if

    call mpi_bcast(kv_top_pressure,   1,   mpi_real8, masterprocid, mpicom, ierr)
    if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
    call mpi_bcast(kv_top_scale,      1,   mpi_real8, masterprocid, mpicom, ierr)
    if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
    call mpi_bcast(kv_freetrop_scale, 1,   mpi_real8, masterprocid, mpicom, ierr)
    if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")

    call mpi_bcast(eddy_lbulk_max,    1,   mpi_real8, masterprocid, mpicom, ierr)
    if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
    call mpi_bcast(eddy_leng_max,     1,   mpi_real8, masterprocid, mpicom, ierr)
    if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
    call mpi_bcast(eddy_max_bot_pressure,     1,   mpi_real8, masterprocid, mpicom, ierr)
    if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")
    call mpi_bcast(eddy_moist_entrain_a2l,    1,   mpi_real8, masterprocid, mpicom, ierr)
    if (ierr /= 0) call endrun(errMsg(__FILE__, __LINE__)//" mpi_bcast error")

  end subroutine eddy_diff_readnl

  subroutine eddy_diff_register()
  
    use physics_buffer, only: pbuf_add_field, dtype_r8, dtype_i4

    character(16) shallow_scheme

    ! Check for UNICON and add relevant pbuf entries.
    call phys_getopts(shallow_scheme_out=shallow_scheme)

    unicon_is_on = (shallow_scheme == "UNICON")

    if (unicon_is_on) then
      call pbuf_add_field('bprod',    'global', dtype_r8, (/pcols,pverp/), bprod_idx)
      call pbuf_add_field('ipbl',     'global', dtype_i4, (/pcols/),       ipbl_idx)
      call pbuf_add_field('kpblh',    'global', dtype_i4, (/pcols/),       kpblh_idx)
      call pbuf_add_field('wstarPBL', 'global', dtype_r8, (/pcols/),       wstarPBL_idx)
      call pbuf_add_field('tkes',     'global', dtype_r8, (/pcols/),       tkes_idx)
      call pbuf_add_field('went',     'global', dtype_r8, (/pcols/),       went_idx)
    end if

  end subroutine eddy_diff_register

  subroutine eddy_diff_init(pbuf2d, ntop_eddy_in, nbot_eddy_in)

    use error_messages  , only: handle_errmsg
    use cam_history     , only: addfld, add_default, horiz_only
    use constituents    , only: cnst_get_ind
    use ref_pres        , only: pref_mid
    use diffusion_solver, only: new_fieldlist_vdiff, vdiff_select
    use eddy_diff       , only: init_eddy_diff
    use physics_buffer  , only: pbuf_set_field, pbuf_get_index

    type(physics_buffer_desc), pointer :: pbuf2d(:,:) ! Physics buffer
    integer, intent(in) :: ntop_eddy_in ! Top interface level to which eddy vertical diffusivity is applied ( = 1 )
    integer, intent(in) :: nbot_eddy_in ! Bottom interface level to which eddy vertical diffusivity is applied ( = pver )

    character(128) errstring

    real(r8) leng_max(pver)
    integer k
    logical history_amwg

    call eddy_diff_final()

    ntop_eddy = ntop_eddy_in
    nbot_eddy = nbot_eddy_in

    do k = 1, pver
      if (pref_mid(k) <= eddy_max_bot_pressure * 1.0e2_r8) then
        leng_max(k) = eddy_leng_max
      else
        leng_max(k) = 40.0e3_r8
      end if
    end do

    if (masterproc) then
      write(iulog, *) 'init_eddy_diff: nturb=', nturb
      write(iulog, *) 'init_eddy_diff: eddy_leng_max=', eddy_leng_max, ' lbulk_max=', eddy_lbulk_max
      do k = 1, pver
        write(iulog, *) 'init_eddy_diff:', k, pref_mid(k), 'leng_max=', leng_max(k)
      end do
    end if

    call init_eddy_diff(pver, gravit, cpair, rair, zvir, &
         latvap, latice, ntop_eddy, nbot_eddy, karman, &
         eddy_lbulk_max, leng_max, &
         eddy_moist_entrain_a2l, errstring)
    call handle_errmsg(errstring, subname='init_eddy_diff')

    ! Set the square of the mixing lengths.
    allocate(ml2(pver+1))
    ml2(1:ntop_eddy) = 0
    do k = ntop_eddy + 1, nbot_eddy
       ml2(k) = 30.0_r8**2
    end do
    ml2(nbot_eddy+1:pver+1) = 0

    ! Get fieldlists to pass to diffusion solver.
    fieldlist_wet   = new_fieldlist_vdiff(1)
    fieldlist_molec = new_fieldlist_vdiff(1)

    call handle_errmsg(vdiff_select(fieldlist_wet,'s'  ), subname='vdiff_select')
    call handle_errmsg(vdiff_select(fieldlist_wet,'q',1), subname='vdiff_select')
    call handle_errmsg(vdiff_select(fieldlist_wet,'u'  ), subname='vdiff_select')
    call handle_errmsg(vdiff_select(fieldlist_wet,'v'  ), subname='vdiff_select')

    ! Cloud mass constituents
    call cnst_get_ind('CLDLIQ', ixcldliq)
    call cnst_get_ind('CLDICE', ixcldice)

    ! Input pbuf fields
    qrl_idx   = pbuf_get_index('QRL')
    wsedl_idx = pbuf_get_index('WSEDL')

    ! Initialize output pbuf fields
    if (is_first_step() .and. unicon_is_on) then
      call pbuf_set_field(pbuf2d, bprod_idx   , 1.0e-5_r8)
      call pbuf_set_field(pbuf2d, ipbl_idx    , 0        )
      call pbuf_set_field(pbuf2d, kpblh_idx   , 1        )
      call pbuf_set_field(pbuf2d, wstarPBL_idx, 0.0_r8   )
      call pbuf_set_field(pbuf2d, tkes_idx    , 0.0_r8   )
      call pbuf_set_field(pbuf2d, went_idx    , 0.0_r8   )
    end if

    ! Scheme-specific default output.
    call phys_getopts(history_amwg_out=history_amwg)

    call addfld('WGUSTD', horiz_only, 'A', 'm/s', 'Wind gusts from turbulence')
    if (history_amwg) then
      call add_default('WGUSTD  ', 1, ' ')
    end if

    ! ------------------------------------------------------------------- !
    ! Writing outputs for detailed analysis of UW moist turbulence scheme !
    ! ------------------------------------------------------------------- !

    call addfld('BPROD'        , [ 'ilev' ], 'A',        'm2/s3', 'Buoyancy Production'     )
    call addfld('SFI'          , [ 'ilev' ], 'A',            '1', 'Interface-layer sat frac')
    call addfld('SPROD'        , [ 'ilev' ], 'A',        'm2/s3', 'Shear Production'        )
    call addfld('UW_errorPBL'  , horiz_only, 'A',         'm2/s', 'Error function of UW PBL')
    call addfld('UW_n2'        , [ 'lev'  ], 'A',          's-2', 'Buoyancy Frequency, LI')
    call addfld('UW_s2'        , [ 'lev'  ], 'A',          's-2', 'Shear Frequency, LI')
    call addfld('UW_ri'        , [ 'lev'  ], 'A',            '1', 'Interface Richardson Number, I')
    call addfld('UW_sfuh'      , [ 'lev'  ], 'A',            '1', 'Upper-Half Saturation Fraction, L')
    call addfld('UW_sflh'      , [ 'lev'  ], 'A',            '1', 'Lower-Half Saturation Fraction, L')
    call addfld('UW_sfi'       , [ 'ilev' ], 'A',            '1', 'Interface Saturation Fraction, I')
    call addfld('UW_cldn'      , [ 'lev'  ], 'A',            '1', 'Cloud Fraction, L')
    call addfld('UW_qrl'       , [ 'lev'  ], 'A', 'gravity W/m2', 'LW cooling rate, L')
    call addfld('UW_ql'        , [ 'lev'  ], 'A',        'kg/kg', 'ql(LWC), L')
    call addfld('UW_chu'       , [ 'ilev' ], 'A', 'gravity kg/J', 'Buoyancy Coefficient, chu, I')
    call addfld('UW_chs'       , [ 'ilev' ], 'A', 'gravity kg/J', 'Buoyancy Coefficient, chs, I')
    call addfld('UW_cmu'       , [ 'ilev' ], 'A','gravity/kg/kg', 'Buoyancy Coefficient, cmu, I')
    call addfld('UW_cms'       , [ 'ilev' ], 'A','gravity/kg/kg', 'Buoyancy Coefficient, cms, I')
    call addfld('UW_tke'       , [ 'ilev' ], 'A',        'm2/s2', 'TKE, I')
    call addfld('UW_wcap'      , [ 'ilev' ], 'A',        'm2/s2', 'Wcap, I')
    call addfld('UW_bprod'     , [ 'ilev' ], 'A',        'm2/s3', 'Buoyancy production, I')
    call addfld('UW_sprod'     , [ 'ilev' ], 'A',        'm2/s3', 'Shear production, I')
    call addfld('UW_kvh'       , [ 'ilev' ], 'A',         'm2/s', 'Eddy diffusivity of heat, I')
    call addfld('UW_kvm'       , [ 'ilev' ], 'A',         'm2/s', 'Eddy diffusivity of uv, I')
    call addfld('UW_pblh'      , horiz_only, 'A',            'm', 'PBLH, 1')
    call addfld('UW_pblhp'     , horiz_only, 'A',           'Pa', 'PBLH pressure, 1')
    call addfld('UW_tpert'     , horiz_only, 'A',            'K', 'Convective T excess, 1')
    call addfld('UW_qpert'     , horiz_only, 'A',        'kg/kg', 'Convective qt excess, I')
    call addfld('UW_wpert'     , horiz_only, 'A',          'm/s', 'Convective W excess, I')
    call addfld('UW_ustar'     , horiz_only, 'A',          'm/s', 'Surface Frictional Velocity, 1')
    call addfld('UW_tkes'      , horiz_only, 'A',        'm2/s2', 'Surface TKE, 1')
    call addfld('UW_minpblh'   , horiz_only, 'A',            'm',  'Minimum PBLH, 1')
    call addfld('UW_turbtype'  , [ 'ilev' ], 'A',            '1',  'Interface Turbulence Type, I')
    call addfld('UW_kbase_o'   , [ 'lev'  ], 'A',            '1',  'Initial CL Base Exterbal Interface Index, CL')
    call addfld('UW_ktop_o'    , [ 'lev'  ], 'A',            '1',  'Initial Top Exterbal Interface Index, CL')
    call addfld('UW_ncvfin_o'  , horiz_only, 'A',            '1',  'Initial Total Number of CL regimes, CL')
    call addfld('UW_kbase_mg'  , [ 'lev'  ], 'A',            '1',  'kbase after merging, CL')
    call addfld('UW_ktop_mg'   , [ 'lev'  ], 'A',            '1',  'ktop after merging, CL')
    call addfld('UW_ncvfin_mg' , horiz_only, 'A',           '1',  'ncvfin after merging, CL')
    call addfld('UW_kbase_f'   , [ 'lev'  ], 'A',            '1',  'Final kbase with SRCL, CL')
    call addfld('UW_ktop_f'    , [ 'lev'  ], 'A',            '1',  'Final ktop with SRCL, CL')
    call addfld('UW_ncvfin_f'  , horiz_only, 'A',            '1',  'Final ncvfin with SRCL, CL')
    call addfld('UW_wet'       , [ 'lev'  ], 'A',          'm/s',  'Entrainment rate at CL top, CL')
    call addfld('UW_web'       , [ 'lev'  ], 'A',          'm/s',  'Entrainment rate at CL base, CL')
    call addfld('UW_jtbu'      , [ 'lev'  ], 'A',         'm/s2',  'Buoyancy jump across CL top, CL')
    call addfld('UW_jbbu'      , [ 'lev'  ], 'A',         'm/s2',  'Buoyancy jump across CL base, CL')
    call addfld('UW_evhc'      , [ 'lev'  ], 'A',            '1',  'Evaporative enhancement factor, CL')
    call addfld('UW_jt2slv'    , [ 'lev'  ], 'A',         'J/kg',  'slv jump for evhc, CL')
    call addfld('UW_n2ht'      , [ 'lev'  ], 'A',          's-2',  'n2 at just below CL top interface, CL')
    call addfld('UW_n2hb'      , [ 'lev'  ], 'A',          's-2',  'n2 at just above CL base interface')
    call addfld('UW_lwp'       , [ 'lev'  ], 'A',        'kg/m2',  'LWP in the CL top layer, CL')
    call addfld('UW_optdepth'  , [ 'lev'  ], 'A',            '1',  'Optical depth of the CL top layer, CL')
    call addfld('UW_radfrac'   , [ 'lev'  ], 'A',            '1',  'Fraction of radiative cooling confined in the CL top')
    call addfld('UW_radf'      , [ 'lev'  ], 'A',        'm2/s3',  'Buoyancy production at the CL top by radf, I')
    call addfld('UW_wstar'     , [ 'lev'  ], 'A',          'm/s',  'Convective velocity, Wstar, CL')
    call addfld('UW_wstar3fact', [ 'lev'  ], 'A',            '1',  'Enhancement of wstar3 due to entrainment, CL')
    call addfld('UW_ebrk'      , [ 'lev'  ], 'A',        'm2/s2',  'CL-averaged TKE, CL')
    call addfld('UW_wbrk'      , [ 'lev'  ], 'A',        'm2/s2',  'CL-averaged W, CL')
    call addfld('UW_lbrk'      , [ 'lev'  ], 'A',            'm',  'CL internal thickness, CL')
    call addfld('UW_ricl'      , [ 'lev'  ], 'A',            '1',  'CL-averaged Ri, CL')
    call addfld('UW_ghcl'      , [ 'lev'  ], 'A',            '1',  'CL-averaged gh, CL')
    call addfld('UW_shcl'      , [ 'lev'  ], 'A',            '1',  'CL-averaged sh, CL')
    call addfld('UW_smcl'      , [ 'lev'  ], 'A',            '1',  'CL-averaged sm, CL')
    call addfld('UW_gh'        , [ 'ilev' ], 'A',            '1',  'gh at all interfaces, I')
    call addfld('UW_sh'        , [ 'ilev' ], 'A',            '1',  'sh at all interfaces, I')
    call addfld('UW_sm'        , [ 'ilev' ], 'A',            '1',  'sm at all interfaces, I')
    call addfld('UW_ria'       , [ 'ilev' ], 'A',            '1',  'ri at all interfaces, I')
    call addfld('UW_leng'      , [ 'ilev' ], 'A',          'm/s',  'Turbulence length scale, I')
    call addfld('UW_wsed'      , [ 'lev'  ], 'A',          'm/s',  'Sedimentation velocity at CL top, CL')

  end subroutine eddy_diff_init

  subroutine eddy_diff_final()

    if (allocated(ml2)) deallocate(ml2)

  end subroutine eddy_diff_final

  subroutine eddy_diff_tend(state, pbuf, cam_in, &
                            ztodt, p, tint, rhoi, cldn, wstarent, &
                            kvm_in, kvh_in, ksrftms, dragblj,tauresx, tauresy, &
                            rrho, ustar, pblh, kvm, kvh, kvq, cgh, cgs, tpert, qpert, &
                            tke, sprod, sfi, turbtype, sm_aw)

    use physics_types , only: physics_state
    use camsrfexch    , only: cam_in_t
    use coords_1d     , only: Coords1D

    type(physics_state)               , intent(in) :: state
    type(physics_buffer_desc), pointer, intent(in) :: pbuf(:)
    type(cam_in_t)                    , intent(in) :: cam_in
    real(r8)                          , intent(in) :: ztodt
    type(Coords1D)                    , intent(in) :: p
    real(r8)                          , intent(in) :: tint(pcols,pver+1)
    real(r8)                          , intent(in) :: rhoi(pcols,pver+1)
    real(r8)                          , intent(in) :: cldn(pcols,pver)
    logical                           , intent(in) :: wstarent
    real(r8)                          , intent(in) :: kvm_in(pcols,pver+1)
    real(r8)                          , intent(in) :: kvh_in(pcols,pver+1)
    real(r8)                          , intent(in) :: ksrftms(pcols)
    real(r8)                          , intent(in) :: dragblj(pcols,pver)       ! Drag profile from Beljaars SGO form drag [ 1/s ]
    real(r8)                          , intent(inout) :: tauresx(pcols)
    real(r8)                          , intent(inout) :: tauresy(pcols)
    real(r8)                          , intent(out) :: rrho(pcols)
    real(r8)                          , intent(out) :: ustar(pcols)
    real(r8)                          , intent(out) :: pblh(pcols)
    real(r8)                          , intent(out) :: kvm(pcols,pver+1)
    real(r8)                          , intent(out) :: kvh(pcols,pver+1)
    real(r8)                          , intent(out) :: kvq(pcols,pver+1)
    real(r8)                          , intent(out) :: cgh(pcols,pver+1)
    real(r8)                          , intent(out) :: cgs(pcols,pver+1)
    real(r8)                          , intent(out) :: tpert(pcols)
    real(r8)                          , intent(out) :: qpert(pcols)
    real(r8)                          , intent(out) :: tke(pcols,pver+1)
    real(r8)                          , intent(out) :: sprod(pcols,pver+1)
    real(r8)                          , intent(out) :: sfi(pcols,pver+1)
    integer(i4)                       , intent(out) :: turbtype(pcols,pver+1)
    real(r8)                          , intent(out) :: sm_aw(pcols,pver+1)

    integer i, k

    call compute_eddy_diff( pbuf, state%lchnk    ,                                     &
         pcols    , pver        , state%ncol       , state%t    , tint, state%q(:,:,1) , ztodt   , &
         state%q(:,:,ixcldliq)  , state%q(:,:,ixcldice)   ,                            &
         state%s  , p           , rhoi, cldn       , &
         state%zm , state%zi    , state%pmid , state%pint , state%u        , state%v , &
         cam_in%wsx, cam_in%wsy , cam_in%shf , cam_in%cflx(:,1) , wstarent           , &
         rrho     , ustar       , pblh       , kvm_in     , kvh_in         , kvm     , &
         kvh      , kvq         , cgh        ,                                         &
         cgs      , tpert       , qpert      , tke            ,                        &
         sprod    , sfi         ,                                                      &
         tauresx  , tauresy     , ksrftms    , dragblj , turbtype   , sm_aw )

    ! The diffusivities from diag_TKE can be much larger than from HB in the free
    ! troposphere and upper atmosphere. These seem to be larger than observations,
    ! and in WACCM the gw_drag code is already applying an eddy diffusivity in the
    ! upper atmosphere. Optionally, adjust the diffusivities in the free troposphere
    ! or the upper atmosphere.
    !
    ! NOTE: Further investigation should be done as to why the diffusivities are
    ! larger in diag_TKE.
    if (kv_freetrop_scale /= 1 .or. (kv_top_scale /= 1 .and. kv_top_pressure > 0)) then
      do i = 1, state%ncol
        do k = 1, pverp
          ! Outside of the boundary layer?
          if (state%zi(i,k) > pblh(i)) then
            ! In the upper atmosphere?
            if (state%pint(i,k) <= kv_top_pressure) then
              kvh(i,k) = kvh(i,k) * kv_top_scale
              kvm(i,k) = kvm(i,k) * kv_top_scale
              kvq(i,k) = kvq(i,k) * kv_top_scale
            else
              kvh(i,k) = kvh(i,k) * kv_freetrop_scale
              kvm(i,k) = kvm(i,k) * kv_freetrop_scale
              kvq(i,k) = kvq(i,k) * kv_freetrop_scale
            end if
          else
            exit
          end if
        end do
      end do
    end if

  end subroutine eddy_diff_tend

  subroutine compute_eddy_diff( pbuf, lchnk  ,                                                      &
                                pcols  , pver   , ncol     , t       , tint, qv       , ztodt   ,   &
                                ql     , qi     , s        , p       , rhoi, cldn     ,             &
                                z      , zi     , pmid     , pi      , u        , v       ,         &
                                taux   , tauy   , shflx    , qflx    , wstarent ,           rrho  , &
                                ustar  , pblh   , kvm_in   , kvh_in  , kvm_out  , kvh_out , kvq   , &
                                cgh    , cgs    , tpert    , qpert   , tke     ,                    &
                                sprod  , sfi    ,                                                   &
                                tauresx, tauresy, ksrftms, dragblj, turbtype, sm_aw )

    !-------------------------------------------------------------------- !
    ! Purpose: Interface to compute eddy diffusivities.                   !
    !          Eddy diffusivities are calculated in a fully implicit way  !
    !          through iteration process.                                 !
    ! Author:  Sungsu Park. August. 2006.                                 !
    !                       May.    2008.                                 !
    !-------------------------------------------------------------------- !

    use diffusion_solver, only: compute_vdiff
    use cam_history     , only: outfld
    use phys_debug_util , only: phys_debug_col
    use physconst       , only: cpairv
    use pbl_utils       , only: calc_ustar, austausch_atm
    use error_messages  , only: handle_errmsg
    use coords_1d       , only: Coords1D
    use wv_saturation   , only: qsat
    use eddy_diff       , only: trbintd, caleddy
    use physics_buffer  , only: pbuf_get_field

    type(physics_buffer_desc), pointer, intent(in) :: pbuf(:)
    integer,  intent(in)    :: lchnk
    integer,  intent(in)    :: pcols                     ! Number of atmospheric columns [ # ]
    integer,  intent(in)    :: pver                      ! Number of atmospheric layers  [ # ]
    integer,  intent(in)    :: ncol                      ! Number of atmospheric columns [ # ]
    logical,  intent(in)    :: wstarent                  ! .true. means use the 'wstar' entrainment closure.
    real(r8), intent(in)    :: ztodt                     ! Physics integration time step 2 delta-t [ s ]
    real(r8), intent(in)    :: t(pcols,pver)             ! Temperature [ K ]
    real(r8), intent(in)    :: tint(pcols,pver+1)        ! Temperature defined on interfaces [ K ]
    real(r8), intent(in)    :: qv(pcols,pver)            ! Water vapor  specific humidity [ kg/kg ]
    real(r8), intent(in)    :: ql(pcols,pver)            ! Liquid water specific humidity [ kg/kg ]
    real(r8), intent(in)    :: qi(pcols,pver)            ! Ice specific humidity [ kg/kg ]
    real(r8), intent(in)    :: s(pcols,pver)             ! Dry static energy [ J/kg ]
    type(Coords1D), intent(in) :: p                      ! Pressure coordinates for solver [ Pa ]
    real(r8), intent(in)    :: rhoi(pcols,pver+1)        ! Density at interfaces [ kg/m3 ]
    real(r8), intent(in)    :: cldn(pcols,pver)          ! Stratiform cloud fraction [ fraction ]
    real(r8), intent(in)    :: z(pcols,pver)             ! Layer mid-point height above surface [ m ]
    real(r8), intent(in)    :: zi(pcols,pver+1)          ! Interface height above surface [ m ]
    real(r8), intent(in)    :: pmid(pcols,pver)          ! Layer mid-point pressure [ Pa ]
    real(r8), intent(in)    :: pi(pcols,pver+1)          ! Interface pressure [ Pa ]
    real(r8), intent(in)    :: u(pcols,pver)             ! Zonal velocity [ m/s ]
    real(r8), intent(in)    :: v(pcols,pver)             ! Meridional velocity [ m/s ]
    real(r8), intent(in)    :: taux(pcols)               ! Zonal wind stress at surface [ N/m2 ]
    real(r8), intent(in)    :: tauy(pcols)               ! Meridional wind stress at surface [ N/m2 ]
    real(r8), intent(in)    :: shflx(pcols)              ! Sensible heat flux at surface [ unit ? ]
    real(r8), intent(in)    :: qflx(pcols)               ! Water vapor flux at surface [ unit ? ]
    real(r8), intent(in)    :: kvm_in(pcols,pver+1)      ! kvm saved from last timestep [ m2/s ]
    real(r8), intent(in)    :: kvh_in(pcols,pver+1)      ! kvh saved from last timestep [ m2/s ]
    real(r8), intent(in)    :: ksrftms(pcols)            ! Surface drag coefficient of turbulent mountain stress [ unit ? ]
    real(r8), intent(in)    :: dragblj(pcols,pver)       ! Drag profile from Beljaars SGO form drag [ 1/s ]

    real(r8), intent(out)   :: kvm_out(pcols,pver+1)     ! Eddy diffusivity for momentum [ m2/s ]
    real(r8), intent(out)   :: kvh_out(pcols,pver+1)     ! Eddy diffusivity for heat [ m2/s ]
    real(r8), intent(out)   :: kvq(pcols,pver+1)         ! Eddy diffusivity for constituents, moisture and tracers [ m2/s ]
                                                         ! (note not having '_out')
    real(r8), intent(out)   :: rrho(pcols)               ! Reciprocal of density at the lowest layer
    real(r8), intent(out)   :: ustar(pcols)              ! Surface friction velocity [ m/s ]
    real(r8), intent(out)   :: pblh(pcols)               ! PBL top height [ m ]
    real(r8), intent(out)   :: cgh(pcols,pver+1)         ! Counter-gradient term for heat [ J/kg/m ]
    real(r8), intent(out)   :: cgs(pcols,pver+1)         ! Counter-gradient star [ cg/flux ]
    real(r8), intent(out)   :: tpert(pcols)              ! Convective temperature excess [ K ]
    real(r8), intent(out)   :: qpert(pcols)              ! Convective humidity excess [ kg/kg ]
    real(r8), intent(out)   :: tke(pcols,pver+1)         ! Turbulent kinetic energy [ m2/s2 ]
    real(r8), intent(out)   :: sprod(pcols,pver+1)       ! Shear production [ m2/s3 ]
    real(r8), intent(out)   :: sfi(pcols,pver+1)         ! Interfacial layer saturation fraction [ fraction ]
    integer(i4), intent(out):: turbtype(pcols,pver+1)    ! Turbulence type identifier at all interfaces [ no unit ]
    real(r8), intent(out)   :: sm_aw(pcols,pver+1)       ! Normalized Galperin instability function for momentum [ no unit ]
                                                         ! This is 1 when neutral condition (Ri=0),
                                                         ! 4.964 for maximum unstable case, and 0 when Ri > Ricrit=0.19.

    real(r8), intent(inout) :: tauresx(pcols)            ! Residual stress to be added in vdiff to correct for turb
    real(r8), intent(inout) :: tauresy(pcols)            ! Stress mismatch between sfc and atm accumulated in prior timesteps

    real(r8), pointer :: qrl(:,:)                        ! LW radiative cooling rate
    real(r8), pointer :: wsedl(:,:)                      ! Sedimentation velocity
                                                         ! of stratiform liquid cloud droplet [ m/s ]

    real(r8), pointer :: bprod(:,:)                      ! Buoyancy production of tke [ m2/s3 ]
    integer(i4), pointer :: ipbl(:)                      ! If 1, PBL is CL, while if 0, PBL is STL.
    integer(i4), pointer :: kpblh(:)                     ! Layer index containing PBL top within or at the base interface
    real(r8), pointer :: wstarPBL(:)                     ! Convective velocity within PBL [ m/s ]
    real(r8), pointer :: tkes(:)                         ! TKE at surface interface [ m2/s2 ]
    real(r8), pointer :: went(:)                         ! Entrainment rate at the PBL top interface [ m/s ]

    integer icol, i, k, iturb, status

    character(2048)         :: warnstring                ! Warning(s) to print
    character(128)          :: errstring                 ! Error message

    real(r8)                :: kvf(pcols,pver+1)         ! Free atmospheric eddy diffusivity [ m2/s ]
    real(r8)                :: kvm(pcols,pver+1)         ! Eddy diffusivity for momentum [ m2/s ]
    real(r8)                :: kvh(pcols,pver+1)         ! Eddy diffusivity for heat [ m2/s ]
    real(r8)                :: kvm_preo(pcols,pver+1)    ! Eddy diffusivity for momentum [ m2/s ]
    real(r8)                :: kvh_preo(pcols,pver+1)    ! Eddy diffusivity for heat [ m2/s ]
    real(r8)                :: kvm_pre(pcols,pver+1)     ! Eddy diffusivity for momentum [ m2/s ]
    real(r8)                :: kvh_pre(pcols,pver+1)     ! Eddy diffusivity for heat [ m2/s ]
    real(r8)                :: errorPBL(pcols)           ! Error function showing whether PBL produced convergent solution or not.
                                                         ! [ unit ? ]
    real(r8)                :: s2(pcols,pver)            ! Shear squared, defined at interfaces except surface [ s-2 ]
    real(r8)                :: n2(pcols,pver)            ! Buoyancy frequency, defined at interfaces except surface [ s-2 ]
    real(r8)                :: ri(pcols,pver)            ! Richardson number, 'n2/s2', defined at interfaces except surface [ s-2 ]
    real(r8)                :: pblhp(pcols)              ! PBL top pressure [ Pa ]
    real(r8)                :: minpblh(pcols)            ! Minimum PBL height based on surface stress

    real(r8)                :: qt(pcols,pver)            ! Total specific humidity [ kg/kg ]
    real(r8)                :: sfuh(pcols,pver)          ! Saturation fraction in upper half-layer [ fraction ]
    real(r8)                :: sflh(pcols,pver)          ! Saturation fraction in lower half-layer [ fraction ]
    real(r8)                :: sl(pcols,pver)            ! Liquid water static energy [ J/kg ]
    real(r8)                :: slv(pcols,pver)           ! Liquid water virtual static energy [ J/kg ]
    real(r8)                :: slslope(pcols,pver)       ! Slope of 'sl' in each layer
    real(r8)                :: qtslope(pcols,pver)       ! Slope of 'qt' in each layer
    real(r8)                :: qvfd(pcols,pver)          ! Specific humidity for diffusion [ kg/kg ]
    real(r8)                :: tfd(pcols,pver)           ! Temperature for diffusion [ K ]
    real(r8)                :: slfd(pcols,pver)          ! Liquid static energy [ J/kg ]
    real(r8)                :: qtfd(pcols,pver)          ! Total specific humidity [ kg/kg ]
    real(r8)                :: qlfd(pcols,pver)          ! Liquid water specific humidity for diffusion [ kg/kg ]
    real(r8)                :: ufd(pcols,pver)           ! U-wind for diffusion [ m/s ]
    real(r8)                :: vfd(pcols,pver)           ! V-wind for diffusion [ m/s ]

    ! Buoyancy coefficients : w'b' = ch * w'sl' + cm * w'qt'

    real(r8)                :: chu(pcols,pver+1)         ! Heat buoyancy coef for dry states, defined at each interface, finally.
    real(r8)                :: chs(pcols,pver+1)         ! Heat buoyancy coef for sat states, defined at each interface, finally.
    real(r8)                :: cmu(pcols,pver+1)         ! Moisture buoyancy coef for dry states,
                                                         ! defined at each interface, finally.
    real(r8)                :: cms(pcols,pver+1)         ! Moisture buoyancy coef for sat states,
                                                         ! defined at each interface, finally.

    real(r8)                :: jnk1d(pcols)
    real(r8)                :: jnk2d(pcols,pver+1)
    real(r8)                :: zero(pcols)
    real(r8)                :: zero2d(pcols,pver+1)
    real(r8)                :: es                     ! Saturation vapor pressure
    real(r8)                :: qs                     ! Saturation specific humidity
    real(r8)                :: ep2, templ, temps

    ! ------------------------------- !
    ! Variables for diagnostic output !
    ! ------------------------------- !

    real(r8)                :: wpert(pcols)              ! Turbulent velocity excess [ m/s ]

    real(r8)                :: kbase_o(pcols,ncvmax)     ! Original external base interface index of CL from 'exacol'
    real(r8)                :: ktop_o(pcols,ncvmax)      ! Original external top  interface index of CL from 'exacol'
    real(r8)                :: ncvfin_o(pcols)           ! Original number of CLs from 'exacol'
    real(r8)                :: kbase_mg(pcols,ncvmax)    ! 'kbase' after extending-merging from 'zisocl'
    real(r8)                :: ktop_mg(pcols,ncvmax)     ! 'ktop' after extending-merging from 'zisocl'
    real(r8)                :: ncvfin_mg(pcols)          ! 'ncvfin' after extending-merging from 'zisocl'
    real(r8)                :: kbase_f(pcols,ncvmax)     ! Final 'kbase' after extending-merging & including SRCL
    real(r8)                :: ktop_f(pcols,ncvmax)      ! Final 'ktop' after extending-merging & including SRCL
    real(r8)                :: ncvfin_f(pcols)           ! Final 'ncvfin' after extending-merging & including SRCL
    real(r8)                :: wet(pcols,ncvmax)         ! Entrainment rate at the CL top  [ m/s ]
    real(r8)                :: web(pcols,ncvmax)         ! Entrainment rate at the CL base [ m/s ].
                                                         ! Set to zero if CL is based at surface.
    real(r8)                :: jtbu(pcols,ncvmax)        ! Buoyancy jump across the CL top  [ m/s2 ]
    real(r8)                :: jbbu(pcols,ncvmax)        ! Buoyancy jump across the CL base [ m/s2 ]
    real(r8)                :: evhc(pcols,ncvmax)        ! Evaporative enhancement factor at the CL top
    real(r8)                :: jt2slv(pcols,ncvmax)      ! Jump of slv ( across two layers ) at CL top used only for evhc [ J/kg ]
    real(r8)                :: n2ht(pcols,ncvmax)        ! n2 defined at the CL top  interface but using
                                                         ! sfuh(kt)   instead of sfi(kt) [ s-2 ]
    real(r8)                :: n2hb(pcols,ncvmax)        ! n2 defined at the CL base interface but using
                                                         ! sflh(kb-1) instead of sfi(kb) [ s-2 ]
    real(r8)                :: lwp(pcols,ncvmax)         ! LWP in the CL top layer [ kg/m2 ]
    real(r8)                :: opt_depth(pcols,ncvmax)   ! Optical depth of the CL top layer
    real(r8)                :: radinvfrac(pcols,ncvmax)  ! Fraction of radiative cooling confined in the top portion of CL top layer
    real(r8)                :: radf(pcols,ncvmax)        ! Buoyancy production at the CL top due to LW radiative cooling [ m2/s3 ]
    real(r8)                :: wstar(pcols,ncvmax)       ! Convective velocity in each CL [ m/s ]
    real(r8)                :: wstar3fact(pcols,ncvmax)  ! Enhancement of 'wstar3' due to entrainment (inverse) [ no unit ]
    real(r8)                :: ebrk(pcols,ncvmax)        ! Net mean TKE of CL including entrainment effect [ m2/s2 ]
    real(r8)                :: wbrk(pcols,ncvmax)        ! Net mean normalized TKE (W) of CL,
                                                         ! 'ebrk/b1' including entrainment effect [ m2/s2 ]
    real(r8)                :: lbrk(pcols,ncvmax)        ! Energetic internal thickness of CL [m]
    real(r8)                :: ricl(pcols,ncvmax)        ! CL internal mean Richardson number
    real(r8)                :: ghcl(pcols,ncvmax)        ! Half of normalized buoyancy production of CL
    real(r8)                :: shcl(pcols,ncvmax)        ! Galperin instability function of heat-moisture of CL
    real(r8)                :: smcl(pcols,ncvmax)        ! Galperin instability function of mementum of CL
    real(r8)                :: ghi(pcols,pver+1)         ! Half of normalized buoyancy production at all interfaces
    real(r8)                :: shi(pcols,pver+1)         ! Galperin instability function of heat-moisture at all interfaces
    real(r8)                :: smi(pcols,pver+1)         ! Galperin instability function of heat-moisture at all interfaces
    real(r8)                :: rii(pcols,pver+1)         ! Interfacial Richardson number defined at all interfaces
    real(r8)                :: lengi(pcols,pver+1)       ! Turbulence length scale at all interfaces [ m ]
    real(r8)                :: wcap(pcols,pver+1)        ! Normalized TKE at all interfaces [ m2/s2 ]
    ! For sedimentation-entrainment feedback
    real(r8)                :: wsed(pcols,ncvmax)        ! Sedimentation velocity at the top of each CL [ m/s ]

    logical,          parameter :: use_kvf        =  .false.      ! .true. (.false.) : initialize kvh/kvm =  kvf ( 0. )
    real(r8),         parameter :: lambda         =   0.5_r8      ! Under-relaxation factor ( 0 < lambda =< 1 )

    zero   = 0
    zero2d = 0

    ! ---------------------------------------------- !
    ! Get LW radiative heating out of physics buffer !
    ! ---------------------------------------------- !
    call pbuf_get_field(pbuf, qrl_idx,   qrl)
    call pbuf_get_field(pbuf, wsedl_idx, wsedl)

    ! These fields are put into the pbuf for UNICON only.
    if (unicon_is_on) then
      call pbuf_get_field(pbuf, bprod_idx,    bprod)
      call pbuf_get_field(pbuf, ipbl_idx,     ipbl)
      call pbuf_get_field(pbuf, kpblh_idx,    kpblh)
      call pbuf_get_field(pbuf, wstarPBL_idx, wstarPBL)
      call pbuf_get_field(pbuf, tkes_idx,     tkes)
      call pbuf_get_field(pbuf, went_idx,     went)
    else
      allocate(bprod(pcols,pverp), ipbl(pcols), kpblh(pcols), wstarPBL(pcols), tkes(pcols), went(pcols))
    end if

    ! ----------------------- !
    ! Main Computation Begins !
    ! ----------------------- !

    ufd (:ncol,:) = u (:ncol,:)
    vfd (:ncol,:) = v (:ncol,:)
    tfd (:ncol,:) = t (:ncol,:)
    qvfd(:ncol,:) = qv(:ncol,:)
    qlfd(:ncol,:) = ql(:ncol,:)

    do iturb = 1, nturb
      ! Total stress includes 'tms'.
      ! Here, in computing 'tms', we can use either iteratively changed 'ufd,vfd' or the
      ! initially given 'u,v' to the PBL scheme. Note that normal stress, 'taux, tauy'
      ! are not changed by iteration. In order to treat 'tms' in a fully implicit way,
      ! I am using updated wind, here.

      ! Compute ustar
      call calc_ustar(ncol, tfd(:ncol,pver), pmid(:ncol,pver), &
                      taux(:ncol) - ksrftms(:ncol) * ufd(:ncol,pver), & ! Zonal wind stress
                      tauy(:ncol) - ksrftms(:ncol) * vfd(:ncol,pver), & ! Meridional wind stress
                      rrho(:ncol), ustar(:ncol))
      minpblh(:ncol) = 100.0_r8 * ustar(:ncol)   ! By construction, 'minpblh' is larger than 1 [m] when 'ustar_min = 0.01'.

      ! Calculate (qt,sl,n2,s2,ri) from a given set of (t,qv,ql,qi,u,v)

      call trbintd(pcols    , pver    , ncol  , z       , ufd     , vfd     , tfd   , pmid    , &
                   s2       , n2      , ri    , zi      , pi      , cldn    , qtfd  , qvfd    , &
                   qlfd     , qi      , sfi   , sfuh    , sflh    , slfd    , slv   , slslope , &
                   qtslope  , chs     , chu   , cms     , cmu     )

      ! Save initial (i.e., before iterative diffusion) profile of (qt,sl) at each iteration.
      ! Only necessary for (qt,sl) not (u,v) because (qt,sl) are newly calculated variables.

      if (iturb == 1) then
        qt(:ncol,:) = qtfd(:ncol,:)
        sl(:ncol,:) = slfd(:ncol,:)
      end if

      ! Get free atmosphere exchange coefficients. This 'kvf' is not used in UW moist PBL scheme
      if (use_kvf) then
        call austausch_atm(pcols, ncol, pver, ntop_eddy, nbot_eddy, ml2, ri, s2, kvf)
      else
        kvf = 0
      end if

      ! Initialize kvh/kvm to send to caleddy, depending on model timestep and iteration number
      ! This is necessary for 'wstar-based' entrainment closure.
      if (iturb == 1) then
        if (is_first_step()) then
          ! First iteration of first model timestep : Use free tropospheric value or zero.
          kvh(:ncol,:) = kvf(:ncol,:)
          kvm(:ncol,:) = kvf(:ncol,:)
        else
          ! First iteration on any model timestep except the first : Use value from previous timestep
          kvh(:ncol,:) = kvh_in(:ncol,:)
          kvm(:ncol,:) = kvm_in(:ncol,:)
        end if
      else
        ! Not the first iteration : Use from previous iteration
        kvh(:ncol,:) = kvh_out(:ncol,:)
        kvm(:ncol,:) = kvm_out(:ncol,:)
      end if

      ! Calculate eddy diffusivity (kvh_out,kvm_out) and (tke,bprod,sprod) using
      ! a given (kvh,kvm) which are used only for initializing (bprod,sprod)  at
      ! the first part of caleddy. (bprod,sprod) are fully updated at the end of
      ! caleddy after calculating (kvh_out,kvm_out)

      call caleddy(pcols     , pver      , ncol      ,                     &
                  slfd      , qtfd      , qlfd      , slv      ,ufd     , &
                  vfd       , pi        , z         , zi       ,          &
                  qflx      , shflx     , slslope   , qtslope  ,          &
                  chu       , chs       , cmu       , cms      ,sfuh    , &
                  sflh      , n2        , s2        , ri       ,rrho    , &
                  pblh      , ustar     ,                                 &
                  kvh       , kvm       , kvh_out   , kvm_out  ,          &
                  tpert     , qpert     , qrl       , kvf      , tke    , &
                  wstarent  , bprod     , sprod     , minpblh  , wpert  , &
                  tkes      , went      , turbtype  , sm_aw    ,          &
                  kbase_o   , ktop_o    , ncvfin_o  ,                     &
                  kbase_mg  , ktop_mg   , ncvfin_mg ,                     &
                  kbase_f   , ktop_f    , ncvfin_f  ,                     &
                  wet       , web       , jtbu      , jbbu     ,          &
                  evhc      , jt2slv    , n2ht      , n2hb     ,          &
                  lwp       , opt_depth , radinvfrac, radf     ,          &
                  wstar     , wstar3fact,                                 &
                  ebrk      , wbrk      , lbrk      , ricl     , ghcl   , &
                  shcl      , smcl      , ghi       , shi      , smi    , &
                  rii       , lengi     , wcap      , pblhp    , cldn   , &
                  ipbl      , kpblh     , wsedl     , wsed, &
                  warnstring, errstring)

      if (trim(warnstring) /= '') then
        write(iulog, *) 'eddy_diff_cam: Messages from caleddy follow.'
        write(iulog, *) warnstring
      end if

      call handle_errmsg(errstring, subname='caleddy')

      ! Calculate errorPBL to check whether PBL produced convergent solutions or not.
      if (iturb == nturb) then
        do i = 1, ncol
          errorPBL(i) = 0._r8
          do k = 1, pver
            errorPBL(i) = errorPBL(i) + ( kvh(i,k) - kvh_out(i,k) )**2
          end do
          errorPBL(i) = sqrt(errorPBL(i)/pver)
        end do
      end if

      ! Eddy diffusivities which will be used for the initialization of (bprod,
      ! sprod) in 'caleddy' at the next iteration step.
      if (iturb > 1 .and. iturb < nturb) then
        kvm_out(:ncol,:) = lambda * kvm_out(:ncol,:) + (1 - lambda) * kvm(:ncol,:)
        kvh_out(:ncol,:) = lambda * kvh_out(:ncol,:) + (1 - lambda) * kvh(:ncol,:)
      end if

      ! Set nonlocal terms to zero for flux diagnostics, since not used by caleddy.
      cgh(:ncol,:) = 0
      cgs(:ncol,:) = 0

      if (iturb < nturb) then
        ! Each time we diffuse the original state

        slfd(:ncol,:)  = sl(:ncol,:)
        qtfd(:ncol,:)  = qt(:ncol,:)
        ufd (:ncol,:)  = u (:ncol,:)
        vfd (:ncol,:)  = v (:ncol,:)

        ! Diffuse initial profile of each time step using a given (kvh_out,kvm_out)
        ! In the below 'compute_vdiff', (slfd,qtfd,ufd,vfd) are 'inout' variables.
        call compute_vdiff(lchnk   ,                                                  &
                            pcols   , pver     , 1        , ncol         , tint, &
                            p       , t        , rhoi     , ztodt        , taux      , &
                            tauy    , shflx    , qflx     , &
                            kvh_out , kvm_out  , kvh_out  , cgs          , cgh       , &
                            zi      , ksrftms  , dragblj  , &
                            zero    , fieldlist_wet, fieldlist_molec, &
                            ufd     , vfd      , qtfd     , slfd         ,             &
                            jnk1d   , jnk1d    , jnk2d    , jnk1d        , errstring , &
                            tauresx , tauresy  , 0        , cpairv(:,:,lchnk), zero, &
                            .false., .false. )

        call handle_errmsg(errstring, subname='compute_vdiff', &
            extra_msg="compute_vdiff called from eddy_diff_cam")

        ! Retrieve (tfd,qvfd,qlfd) from (slfd,qtfd) in order to use 'trbintd' at the next iteration.
        do k = 1, pver
          do i = 1, ncol
            ! ----------------------------------------------------- !
            ! Compute the condensate 'qlfd' in the updated profiles !
            ! ----------------------------------------------------- !
            ! Option.1 : Assume grid-mean condensate is homogeneously diffused by the moist turbulence scheme.
            !            This should be used if 'pseudodiff = .false.' in vertical_diffusion.F90.
            ! Modification : Need to be check whether below is correct in the presence of ice, qi.
            !                I should understand why the variation of ice, qi is neglected during diffusion.
            templ     = (slfd(i,k) - gravit * z(i,k)) / cpair
            call qsat(templ, pmid(i,k), es, qs)
            ep2       = 0.622_r8
            temps     = templ + ( qtfd(i,k) - qs ) / ( cpair / latvap + latvap * qs / ( rair * templ**2 ) )
            call qsat(temps, pmid(i,k), es, qs)
            qlfd(i,k) = max( qtfd(i,k) - qi(i,k) - qs ,0._r8 )
            ! Option.2 : Assume condensate is not diffused by the moist turbulence scheme.
            !            This should bs used if 'pseudodiff = .true.'  in vertical_diffusion.F90.
            ! qlfd(i,k) = ql(i,k)
            ! ----------------------------- !
            ! Compute the other 'qvfd, tfd' !
            ! ----------------------------- !
            qvfd(i,k) = max(0.0_r8, qtfd(i,k) - qi(i,k) - qlfd(i,k))
            tfd(i,k)  = (slfd(i,k) + latvap * qlfd(i,k) + (latvap+latice) * qi(i,k) - gravit*z(i,k)) / cpair
          end do
        end do
      end if
    end do  ! End of 'iturb' iteration

    kvq(:ncol,:) = kvh_out(:ncol,:)

    ! Compute 'wstar' within the PBL for use in the future convection scheme.
    do i = 1, ncol
      if (ipbl(i) == 1) then
        wstarPBL(i) = max(0.0_r8, wstar(i,1))
      else
        wstarPBL(i) = 0
      end if
    end do

    ! --------------------------------------------------------------- !
    ! Writing for detailed diagnostic analysis of UW moist PBL scheme !
    ! --------------------------------------------------------------- !

    call outfld( 'WGUSTD'       , wpert     , pcols, lchnk)
    call outfld( 'BPROD'        , bprod     , pcols, lchnk)
    call outfld( 'SPROD'        , sprod     , pcols, lchnk)
    call outfld( 'SFI'          , sfi       , pcols, lchnk)
    call outfld( 'UW_errorPBL'  , errorPBL  , pcols, lchnk)
    call outfld( 'UW_n2'        , n2        , pcols, lchnk)
    call outfld( 'UW_s2'        , s2        , pcols, lchnk)
    call outfld( 'UW_ri'        , ri        , pcols, lchnk)
    call outfld( 'UW_sfuh'      , sfuh      , pcols, lchnk)
    call outfld( 'UW_sflh'      , sflh      , pcols, lchnk)
    call outfld( 'UW_sfi'       , sfi       , pcols, lchnk)
    call outfld( 'UW_cldn'      , cldn      , pcols, lchnk)
    call outfld( 'UW_qrl'       , qrl       , pcols, lchnk)
    call outfld( 'UW_ql'        , qlfd      , pcols, lchnk)
    call outfld( 'UW_chu'       , chu       , pcols, lchnk)
    call outfld( 'UW_chs'       , chs       , pcols, lchnk)
    call outfld( 'UW_cmu'       , cmu       , pcols, lchnk)
    call outfld( 'UW_cms'       , cms       , pcols, lchnk)
    call outfld( 'UW_tke'       , tke       , pcols, lchnk)
    call outfld( 'UW_wcap'      , wcap      , pcols, lchnk)
    call outfld( 'UW_bprod'     , bprod     , pcols, lchnk)
    call outfld( 'UW_sprod'     , sprod     , pcols, lchnk)
    call outfld( 'UW_kvh'       , kvh_out   , pcols, lchnk)
    call outfld( 'UW_kvm'       , kvm_out   , pcols, lchnk)
    call outfld( 'UW_pblh'      , pblh      , pcols, lchnk)
    call outfld( 'UW_pblhp'     , pblhp     , pcols, lchnk)
    call outfld( 'UW_tpert'     , tpert     , pcols, lchnk)
    call outfld( 'UW_qpert'     , qpert     , pcols, lchnk)
    call outfld( 'UW_wpert'     , wpert     , pcols, lchnk)
    call outfld( 'UW_ustar'     , ustar     , pcols, lchnk)
    call outfld( 'UW_tkes'      , tkes      , pcols, lchnk)
    call outfld( 'UW_minpblh'   , minpblh   , pcols, lchnk)
    call outfld( 'UW_kbase_o'   , kbase_o   , pcols, lchnk)
    call outfld( 'UW_ktop_o'    , ktop_o    , pcols, lchnk)
    call outfld( 'UW_ncvfin_o'  , ncvfin_o  , pcols, lchnk)
    call outfld( 'UW_kbase_mg'  , kbase_mg  , pcols, lchnk)
    call outfld( 'UW_ktop_mg'   , ktop_mg   , pcols, lchnk)
    call outfld( 'UW_ncvfin_mg' , ncvfin_mg , pcols, lchnk)
    call outfld( 'UW_kbase_f'   , kbase_f   , pcols, lchnk)
    call outfld( 'UW_ktop_f'    , ktop_f    , pcols, lchnk)
    call outfld( 'UW_ncvfin_f'  , ncvfin_f  , pcols, lchnk)
    call outfld( 'UW_wet'       , wet       , pcols, lchnk)
    call outfld( 'UW_web'       , web       , pcols, lchnk)
    call outfld( 'UW_jtbu'      , jtbu      , pcols, lchnk)
    call outfld( 'UW_jbbu'      , jbbu      , pcols, lchnk)
    call outfld( 'UW_evhc'      , evhc      , pcols, lchnk)
    call outfld( 'UW_jt2slv'    , jt2slv    , pcols, lchnk)
    call outfld( 'UW_n2ht'      , n2ht      , pcols, lchnk)
    call outfld( 'UW_n2hb'      , n2hb      , pcols, lchnk)
    call outfld( 'UW_lwp'       , lwp       , pcols, lchnk)
    call outfld( 'UW_optdepth'  , opt_depth , pcols, lchnk)
    call outfld( 'UW_radfrac'   , radinvfrac, pcols, lchnk)
    call outfld( 'UW_radf'      , radf      , pcols, lchnk)
    call outfld( 'UW_wstar'     , wstar     , pcols, lchnk)
    call outfld( 'UW_wstar3fact', wstar3fact, pcols, lchnk)
    call outfld( 'UW_ebrk'      , ebrk      , pcols, lchnk)
    call outfld( 'UW_wbrk'      , wbrk      , pcols, lchnk)
    call outfld( 'UW_lbrk'      , lbrk      , pcols, lchnk)
    call outfld( 'UW_ricl'      , ricl      , pcols, lchnk)
    call outfld( 'UW_ghcl'      , ghcl      , pcols, lchnk)
    call outfld( 'UW_shcl'      , shcl      , pcols, lchnk)
    call outfld( 'UW_smcl'      , smcl      , pcols, lchnk)
    call outfld( 'UW_gh'        , ghi       , pcols, lchnk)
    call outfld( 'UW_sh'        , shi       , pcols, lchnk)
    call outfld( 'UW_sm'        , smi       , pcols, lchnk)
    call outfld( 'UW_ria'       , rii       , pcols, lchnk)
    call outfld( 'UW_leng'      , lengi     , pcols, lchnk)
    call outfld( 'UW_wsed'      , wsed      , pcols, lchnk)
    call outfld( 'UW_turbtype'  , real(turbtype,r8), pcols, lchnk )

    if (.not. unicon_is_on) then
      deallocate(bprod, ipbl, kpblh, wstarPBL, tkes, went)
    end if

  end subroutine compute_eddy_diff

end module eddy_diff_cam
